R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

#install.packages("corrplot")
#install.packages("Rfast")
#install.packages("phyloseq")
#install.packages("glmnet")
#install.packages("randcorr")
#install.packages("hrbrthemes")
library(randcorr)
library(ggplot2)
library(hrbrthemes)
library(Rfast)
library(corrplot)

# 1a
create_TC = function(interval, one_time, start_time){
  
  # Create initial matrix
  TC = matrix(0, 240, 1)
  
  # Counts for how many ones have been added / not added
  on_count = 1
  off_count = 0
  
  for(i in 1:240){
    
    # Check if the current time in source is before our start time
    if (i <= start_time) {
      next
    }
    
    # Change value in matrix if it is within the time period to do so 
    # Otherwise increment the time that one has not been added
    # Reset if at the end of a interval
    if(on_count <= one_time) {
      TC[i] = 1
      on_count = on_count + 1
    } else{
      off_count = off_count + 1
      if (off_count == interval - one_time) {
        on_count = 1
        off_count = 0
      }
  
    }
  
  }
  return(TC)
}

# 1a
TC1 = create_TC(30, 15, 0)
TC2 = create_TC(45, 20, 20)
TC3 = create_TC(60, 25, 0)
TC4 = create_TC(40, 15, 0)
TC5 = create_TC(40, 20, 0)
TC6 = create_TC(40, 25, 0)

# 1a
standard_TC1 = scale(TC1)
standard_TC2 = scale(TC2)
standard_TC3 = scale(TC3)
standard_TC4 = scale(TC4)
standard_TC5 = scale(TC5)
standard_TC6 = scale(TC6)
full_TC = matrix(c(standard_TC1, standard_TC2, standard_TC3, standard_TC4, standard_TC5, standard_TC6), 240, 6)

#normalised_TC1 = sqrt(sum(TC1^2))
#normalised_TC2 = sqrt(sum(TC2^2))
#normalised_TC3 = sqrt(sum(TC3^2))
#normalised_TC4 = sqrt(sum(TC4^2))
#normalised_TC5 = sqrt(sum(TC5^2))
#normalised_TC6 = sqrt(sum(TC6^2))

#(TC1 / normalised_TC1)[1:50]
#(TC5 / normalised_TC5)[1:50]
#(TC6 / normalised_TC6)[1:50]

# 1a
par(mfrow=c(2,3))

plot(full_TC[1:nrow(full_TC), 1], ylab = "Temporal source 1", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 2], ylab = "Temporal source 2", pch = 1, type = "line")

plot(full_TC[1:nrow(full_TC), 3], ylab = "Temporal source 3", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 4], ylab = "Temporal source 4", pch = 1, type = "line")

plot(full_TC[1:nrow(full_TC), 5], ylab = "Temporal source 5", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 6], ylab = "Temporal source 6", pch = 1, type = "line")

# If we normalise, the values of TC will be same where as standardised values 
# Are not the same inter TC eg
# All TC will have 0 values still

# 1b

correlation = cor(full_TC)
corrplot.mixed(correlation, tl.srt = 0)

correlation

              [,1]          [,2]          [,3]       [,4]          [,5]       [,6]
[1,]  1.000000e+00  3.175488e-21  1.690309e-01 0.08606630 -4.536411e-22  0.0860663
[2,]  3.175488e-21  1.000000e+00 -2.857143e-02 0.13093073  1.769200e-20 -0.1309307
[3,]  1.690309e-01 -2.857143e-02  1.000000e+00 0.04364358  4.082770e-21  0.1309307
[4,]  8.606630e-02  1.309307e-01  4.364358e-02 1.00000000  7.745967e-01  0.6000000
[5,] -4.536411e-22  1.769200e-20  4.082770e-21 0.77459667  1.000000e+00  0.7745967
[6,]  8.606630e-02 -1.309307e-01  1.309307e-01 0.60000000  7.745967e-01  1.0000000

# Comment on if on TC correlation or not
# 4/5, 5/6 highly correlated as size of circle is size of correlation and relatively positvely correlated

# 1c
create_tmpSM = function(vert_start,vert_end, horo_start, horo_end) {
  tmpSM = matrix(0, 21, 21)
  tmpSM[vert_start:vert_end, horo_start:horo_end] = 1
  return(tmpSM)
}

# 1c
tmpSM1 = create_tmpSM(2, 6, 2, 6)
tmpSM2 = create_tmpSM(2, 6, 15, 19)
tmpSM3 = create_tmpSM(8, 13, 2, 6)
tmpSM4 = create_tmpSM(8, 13, 15, 19)
tmpSM5 = create_tmpSM(15, 19, 2, 6)
tmpSM6 = create_tmpSM(15, 19, 15, 19)

create_df = function(vert_start,vert_end, horo_start, horo_end) {
  df1 = expand.grid(seq_len(21), seq_len(21))
  tmp = matrix(0,21,21)
  tmp[vert_start:vert_end, horo_start:horo_end] = 1
  column = as.vector(tmp)
  df1$data = (column)
  return(df1)
}

plot1 = create_df(2, 6, 2, 6)
plot2 = create_df(2, 6, 15, 19)
plot3 = create_df(8, 13, 2, 6)
plot4 = create_df(8, 13, 15, 19)
plot5 = create_df(15, 19, 2, 6)
plot6 = create_df(15, 19, 15, 19)

ggplot(plot1, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

ggplot(plot2, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

ggplot(plot3, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

ggplot(plot4, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

ggplot(plot5, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

ggplot(plot6, aes(Var1, Var2, fill= data)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

# 1c
tmpSM1 = matrix(tmpSM1, 1, 441)
tmpSM2 = matrix(tmpSM2, 1, 441)
tmpSM3 = matrix(tmpSM3, 1, 441)
tmpSM4 = matrix(tmpSM4, 1, 441)
tmpSM5 = matrix(tmpSM5, 1, 441)
tmpSM6 = matrix(tmpSM6, 1, 441)
tmpSM = matrix(c(tmpSM1, tmpSM2, tmpSM3, tmpSM4, tmpSM5, tmpSM6), 441, 6)
tmpSM = t(tmpSM)

correlation2 = cor(t(tmpSM))
corrplot.mixed(correlation2, upper = "circle")

correlation2

            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
[1,]  1.00000000 -0.06009615 -0.06623127 -0.06623127 -0.06009615 -0.06009615
[2,] -0.06009615  1.00000000 -0.06623127 -0.06623127 -0.06009615 -0.06009615
[3,] -0.06623127 -0.06623127  1.00000000 -0.07299270 -0.06623127 -0.06623127
[4,] -0.06623127 -0.06623127 -0.07299270  1.00000000 -0.06623127 -0.06623127
[5,] -0.06009615 -0.06009615 -0.06623127 -0.06623127  1.00000000 -0.06009615
[6,] -0.06009615 -0.06009615 -0.06623127 -0.06623127 -0.06009615  1.00000000

# Are theres SM all independent? - yes as not correlated
# Why standardisation not important - If we are looking to replicate it then 
# We still just have two main areas - zero and non zero so the scale of non 0
# Doesn't matter as much

#1d
set.seed(30034)
noise_t = rnorm(6 * 240, mean = 0, sd = sqrt(0.25))
noise_s = rnorm(6 * 441, mean = 0, sd = sqrt(0.15))


noise_s = matrix(noise_s, 6, 441)

noise_t = matrix(noise_t, 240, 6)

# 1d
correlation3 = cor(noise_t)
corrplot(correlation3, tl.srt = 0)

correlation4 = cor(t(noise_s))
corrplot(correlation4, tl.srt = 0)

# Check if these two above are correlated
# Highest is like abs(0.12,0.15) so not really that strongly correlated

# 1d
# Both look pretty normally distributed with over 95% of the data covered within -1.96, 1.96
hist(noise_s)

hist(noise_t)

  # Both are equal to 1 even with formula variance = 1.96sd, noise_s = 0.945 so pre close, noise_t slight less 0.936 makes sense with graph
#length(noise_s[(0 - 2*sqrt(0.15) <= noise_s) & (noise_s <= 0 + 2*sqrt(0.15))]) / length(noise_s)
#length(noise_t[(0 - 2*sqrt(0.25) <= noise_t) & (noise_t <= 0 + 2*sqrt(0.25))]) / length(noise_t)

# noise s looks way more normally distributed, while t looks kinda skewed.
# do conf int to normal dist
# Check whether fills var of 1.96*sd (0.95 data)

correlation5 = cor(noise_t %*% noise_s)
#corrplot(correlation5, tl.srt = 0)
subset = correlation5[1:10, 1:10]
corrplot(subset, tl.srt = 0)

subset

             [,1]        [,2]        [,3]        [,4]       [,5]        [,6]       [,7]       [,8]       [,9]       [,10]
 [1,]  1.00000000 -0.08401539  0.08121198 -0.48582752 -0.2122397 -0.58300192  0.7469938  0.1747322 -0.4616947  0.69749810
 [2,] -0.08401539  1.00000000 -0.36170166 -0.08651628  0.3311008 -0.65485963  0.1186072 -0.6562960  0.2556765 -0.27880412
 [3,]  0.08121198 -0.36170166  1.00000000 -0.61666983 -0.2739593  0.05563854  0.4470863  0.6308545 -0.6288561 -0.05763502
 [4,] -0.48582752 -0.08651628 -0.61666983  1.00000000  0.2248447  0.47069260 -0.5268018 -0.4475564  0.4159525 -0.37703467
 [5,] -0.21223971  0.33110081 -0.27395932  0.22484470  1.0000000  0.24773797 -0.1582580 -0.5139208  0.6008665 -0.70374353
 [6,] -0.58300192 -0.65485963  0.05563854  0.47069260  0.2477380  1.00000000 -0.6314982  0.2674007  0.3960874 -0.40361348
 [7,]  0.74699385  0.11860722  0.44708627 -0.52680181 -0.1582580 -0.63149817  1.0000000  0.2406842 -0.5654764  0.25340041
 [8,]  0.17473222 -0.65629601  0.63085453 -0.44755642 -0.5139208  0.26740068  0.2406842  1.0000000 -0.2328050  0.37238708
 [9,] -0.46169469  0.25567646 -0.62885608  0.41595249  0.6008665  0.39608744 -0.5654764 -0.2328050  1.0000000 -0.39341424
[10,]  0.69749810 -0.27880412 -0.05763502 -0.37703467 -0.7037435 -0.40361348  0.2534004  0.3723871 -0.3934142  1.00000000

#corrplot(correlation5, tl.srt = 0)[1:10, 1:10]
# Definitely some correlation eg higher covariance values such as 0.7

#1e

#Can full_tc%*% noise_S and vice versa exist and what happens to them as if we keep we cannot fit model?
X = (full_TC + noise_t) %*% (tmpSM + noise_s)
scaled_X = matrix(NA, 240, 441)
for (i in 1:441) {
  scaled_X[,i] = scale(X[,i])
}


# Plot graph of sample Time series
set.seed(1)
mysample <- sample(nrow(X), 100)
mysample

  [1]  68 167 129 162 215  43  14 210 187  51 225  85  21 106 182  74   7  73  79 213  37 105 217 110 165  34 227 126  89 172 207
 [32]  33  84 163  70 230  42 166 111 148 156  20  44 121  87 206 197  40 211  25 119 122  39 179 204 134  24 160 234 130  45 146
 [63]  22 115 104 195 175 176 103  75  13 193 212  23 218 208 192  29 141 150 200 108  48 209 180 149  31 102 145 223 196  83 118
 [94]  90 158 107  64 164  60 231

sample = X[1:nrow(X), mysample]

#for (j in 1:10) {
#  if (j == 1) {
#    plot(sample[1:nrow(X), (i - 1)*10 + j], type = "line")
#  }
#  else {
#    if (j == 2) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'red')
#    }
#    
#    if (j == 3) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'blue')
#    }
#    
#    if (j == 4) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'green')
#    }
#    
#    if (j == 5) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'yellow')
#    }
#    if (j == 6) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'darkorange')
#    }
#    
#    if (j == 7) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'cyan')
#    }
#    if (j == 8) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'brown')
#    }
#    if (j == 9) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'grey')
#    }
#    
#    if (j == 10) {
#      lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'beige')
#    }
#    
#  }
#  
#}

# What info does the plot of vars give
variances = colVars(X)
plot(variances, type = "line")

NA
NA
NA

# 2a
A_LSR = solve(t(full_TC) %*% full_TC) %*% t(full_TC) %*% scaled_X
D_LSR = scaled_X%*% abs(t(A_LSR))


#  Plot and heatmap together
plot(D_LSR[1:nrow(D_LSR),1], type = "line")

# produce legend for heatmap
#heatmap(matrix(A[1,1:ncol(A_LSR)], 21, 21), ncol, Colv = NA, Rowv = NA, scale="column")



for (num in 1:6) {
  df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
  count = 1
  col = abs(A_LSR[num,])
  for (i in 1:21) {
    for (j in 1:21) {
      df[count, 1] = i
      df[count, 2] = j
      df[count, 3] = col[(i-1) * 21 + j]
      count = count + 1
    }
    
  }
  ggplot(df, aes(V1, V2, fill= V3)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()
}






# WHy linear relo bet 3rd col d_lsr and 30th of X
# Secondly if you look at the slice number 3 build up by third SM. Now read this 3rd SM column-wise, you will find that 30th (2nd column 9th row =30th pixel) pixel position was filled by this 3rd SM and therefore third TC is the only time course that constructs 30th column of X and this is the TC that you retrieve in 3rd Dlsr. There might be some variance in the scatter plot but you will be able to capture visually the linear trend between 3rd Dlsr and 30th X.
#linear relo - 
plot(D_LSR[1:nrow(D_LSR), 3], scaled_X[1:nrow(scaled_X), 30], xlab = "3rd column of D(LSR)", ylab = "30th column of scaled X")

# No linear relo
plot(D_LSR[1:nrow(D_LSR), 4], scaled_X[1:nrow(scaled_X), 30], xlab = "3rd column of D(LSR)", ylab = "30th column of scaled X")

# Gets the sum of the maximum absolute values in the correlation matrix
max_abs_corr_sum = function(mat1, mat2) {
  corre_mat = cor(mat1, mat2)
  c = get_max_correlations(corre_mat)
  return(sum(c))
}

# Returns the values of the maximum absolute values for each row
get_max_correlations = function(cor_mat) {
  c_var = c(rep(0,6))
  abs_mat = abs(cor_mat)
  cols = max.col(abs_mat)
  
  for (i in 1:6) {
    c_var[i] = cor_mat[i, cols[i]]
  }
  return(c_var)
}

#2b
V = 441
MSE = matrix(0, 1, 2)
for (num in seq(from = 0, to = 1, by = 0.01)) {
  lamda = num
  lamda_hat = lamda * V
  A_RR = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * lamda_hat) %*% t(full_TC) %*% scaled_X
  mse = (sum((scaled_X - full_TC %*% abs(A_RR))^2)) + lamda_hat * sum(A_RR^2)
  MSE = rbind(MSE, c(lamda, mse))
}

plot(MSE[2:nrow(MSE), 1], MSE[2:nrow(MSE), 2])

# From the graph it looks like best value of lamda is 1
A_RR = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * 1* V) %*% t(full_TC) %*% scaled_X
D_RR = scaled_X%*% abs(t(A_RR))
c_tlsr = max_abs_corr_sum(full_TC, D_LSR)
c_trr = max_abs_corr_sum(full_TC, D_RR)

# 4.407241
c_tlsr

[1] 4.506323

# 4.47653
c_trr

[1] 4.682376

c_trr - c_tlsr

[1] 0.1760527

lamda = 1000
lamda_hat = lamda * V
A_1000 = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * lamda_hat) %*% t(full_TC) %*% scaled_X


# As A_LSR decrease, A_RR values decrease
plot(abs(A_1000[1:nrow(A_1000), 1]), abs(A_LSR[1:nrow(A_LSR), 1]), xlab = "A_RR", ylab = "A_LSR")

create_noise_matrix = function(seed) {
  set.seed(seed = seed)
  noise_t = rnorm(6 * 240, mean = 0, sd = sqrt(0.25))
  noise_s = rnorm(6 * 441, mean = 0, sd = sqrt(0.15))


  noise_s = matrix(noise_s, 6, 441)
  noise_t = matrix(noise_t, 240, 6)
  
  X = (full_TC + noise_t) %*% (tmpSM + noise_s)
  new_X = scale(X)
  return(new_X)
}
noise_1 = create_noise_matrix(30035)
noise_2 = create_noise_matrix(30037)
noise_3 = create_noise_matrix(30038)
noise_4 = create_noise_matrix(30036)
noise_5 = create_noise_matrix(30039)
noise_6 = create_noise_matrix(30040)
noise_7 = create_noise_matrix(30041)
noise_8 = create_noise_matrix(30042)
noise_9 = create_noise_matrix(30043)
noise_10 = create_noise_matrix(30044)


noise_df = list(noise_1, noise_9, noise_8, noise_7, noise_6, noise_5,
                noise_2, noise_3, noise_4, noise_10
                )
#for(mat in noise_df) {
#  print(mat[1:5, 1:3]
#}

#2c
LR_MSE = matrix(NA, 21, 10)
N = 240
x1 = 21
x2 = 21
nsrcs = 6
V = x1 * x2
i = 1
j = 1
step <- 1/(norm(full_TC %*% t(full_TC)) * 1.1)

for (num in seq(0,1,0.05)) {
  for (x_rep in noise_df) {
    thr <- num*N*step
    Ao <- matrix(0, nsrcs, 1)
    A <- matrix(0, nsrcs, 1)
    Alr <- matrix(0, nsrcs, x1*x2)
    for (k in 1:(x1*x2)) {
      A <- Ao+step*(t(full_TC) %*% (x_rep[,k]-(full_TC%*%Ao)))
      A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
      
      for (rep in 1:10) {
        Ao <- A
        A <- Ao+step * (t(full_TC)%*%(x_rep[,k]-(full_TC%*%Ao)))
        A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
      }
      
      Alr[,k] <- A
    }
    Dlr = x_rep %*% abs(t(Alr))
    LR_MSE[i,j] = sum(sum((x_rep-Dlr%*%abs(Alr))^2)) / (N*V)
    j = j + 1
  }
  j = 1
  i = i + 1
}

LR_avgs = rowMeans(LR_MSE)
plot( seq(0,1, 0.05), LR_avgs, ylab = "MSE(LR)")

# Min value is 0.580 when p = 0.55, at same value it diverges

# Since first index is at 0, we want the 12th value between 0 and 1 which is 0.6 
# So we need to index the 13th value
index = which.min(LR_avgs)
rho = (index - 1) * 0.05

# 2d

thr <- rho*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
Alr <- matrix(0, nsrcs, x1*x2)
for (k in 1:(x1*x2)) {
  A <- Ao+step*(t(full_TC) %*% (scaled_X[,k]-(full_TC%*%Ao)))
  A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  
  for (rep in 1:10) {
    Ao <- A
    A <- Ao+step * (t(full_TC)%*%(scaled_X[,k]-(full_TC%*%Ao)))
    A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
  }
  
  Alr[,k] <- A
}
Dlr = scaled_X %*% abs(t(Alr))

ctrr = max_abs_corr_sum(full_TC, D_RR)
csrr = max_abs_corr_sum(t(tmpSM), t(A_RR))
ctlr = max_abs_corr_sum(full_TC, Dlr)
cslr= max_abs_corr_sum(t(tmpSM), t(Alr))

# Both are true
ctlr - ctrr #= 0.7413487

[1] 0.6743216

cslr - csrr #= 0.6595151

[1] 0.6778317

ctlr #= 5.217879

[1] 5.356698

ctrr #= 4.47653

[1] 4.682376

cslr #= 2.953101

[1] 2.937037

csrr #= 2.293586

[1] 2.259206

plot(D_RR[, 1:4])

plot(Dlr[, 1:4])

plot(Alr[1:4,])

plot(A_RR[1:4, ])

par(mfrow = c(1,2))

for (num in 1:6) {
  plot(D_RR[, 6], type = "line")
  plot(Dlr[, 6], type = "line")
}

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
count = 1
vec = abs(A_RR[num,])
for (num in 1:6) {
  for (i in 1:21) {
    for (j in 1:21) {
      df[count, 1] = i
      df[count, 2] = j
      df[count, 3] = vec[(i-1) * 21 + j]
      count = count + 1
    }
  
  }
}

df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
count = 1
vec = abs(Alr[num,])
for (num in 1:6) {
  for (i in 1:21) {
    for (j in 1:21) {
      df[count, 1] = i
      df[count, 2] = j
      df[count, 3] = vec[(i-1) * 21 + j]
      count = count + 1
    }
  
  }
}

ggplot(df, aes(V1, V2, fill= V3)) + 
  geom_tile() + scale_fill_gradient(low="blue", high="red") +
  theme_ipsum()

par(mfrow=c(1,2))
# 2e
#svd(full_TC)
plot(svd(full_TC)$d, ylab = "Eigenvalue")
#PC 6 smallest

plot(svd(full_TC)$u[1:240, 6])

plot(full_TC[1:240, 6])

plot(svd(full_TC)$u[1:240, 5])

plot(full_TC[1:240, 5])

plot(svd(full_TC)$u[1:240, 4])

plot(full_TC[1:240, 4])

plot(svd(full_TC)$u[1:240, 3])

plot(full_TC[1:240, 3])

plot(svd(full_TC)$u[1:240, 2])

plot(full_TC[1:240, 2])

plot(svd(full_TC)$u[1:240, 1])

plot(full_TC[1:240, 1])

NA
NA